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Abstract 

The genetic regulatory network, which is constructed from 
the time-courses data sets, is always described as highly 
nonlinear differential equations. Mathematical and 
computational modeling technologies focus on efficiently 
identifying the parameters of the nonlinear dynamic 
biological system. Various derivative-free and derivative- 
based optimization technologies have been proposed 
recently to infer the parameters of the S-type genetic 
regulatory networks (S-systems). The S-system is described 
as coupled power-law functions. As the involved genes 
and/or proteins increase, the identification becomes 
increasingly difficult; multiple attractors exist in the system. 
How to develop an optimization algorithm to reduce the 
computation time while keeping the accuracy is necessary. 
In this study, a gradient-based metaheuristics is proposed. 
The computational method starts with the hill-climbing 
optimization, and solves the stagnation phenomenon by 
using a differential climbing operation and migration 
synchronous evolution. This method was tested with four 
biological systems. To show the performance in the solution 
quality and the computation time, we let the learning be 
implemented in a wide search space ([0, 100] for rate 
constants and [-100, 100] for kinetic orders) and initialized all 
parameters at a bad point (the neighbourhood of 80). 
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Introduction 

Data-driven modelling is a corner stone of systems 
biology (Vilela et al., 2008). Researchers have to model 
biological networks such that the regulatory 
interactions between the genes and/or proteins are 
identified. The S-system structure (Savageau, 1976; 
Voit, 2000) is a popular nonlinear model which is able 
to capture the dynamic behaviour of gene regulation 
networks, metabolic pathways or signal transduction 
cascades. This structure directly identifies the 
interaction between genes and/or proteins, and 
possesses good generalization characteristics. 


However, the identification challenges researchers due 
to variables in multimodal distribution. 

Wang et al. (2010) considered two extreme cases to 
determine the parameters' ranges and mean values. 
Xang et al. (2012) fixed the efflux function and 
determined the unknown parameters of the influx 
function through the slope error. Various derivative- 
free-based technologies have recently proposed. Cho 
et al. (2006) used genetic programming. Wang and 
collaborators (Liu and Wang, 2009; Wang and Liu, 
2010) used hybrid differential evolution. Chen et al. 
(2010) hybridized genetic algorithm and simulated 
annealing. Xu et al. (2007) used the neural network 
with the particle swarm optimization. 

Some researchers adopted derivative-based 
technologies to identify the S-systems' parameters. 
Marino and Voit (2006) proposed a step-by-step 
progress in model complexity. Chou et al. (2006) used 
alternative-regression methods. Vilela et al. (2008) 
used eigenvector methods. Kutalik et al. (2007) used 
Newton-flow methods. Sriyudthsak et al. (2013) 
introduced the Granger causality test to infer the 
structure and then used Levenberg-Marquardt 
algorithm to solve the unknown parameters. 
Chemmangattuvalappil et al. (2012) sequentially 
reduced the number of the unknown parameters and 
used least square optimization for parameter 
identification. However, the derivative-based 
optimization methods have the possibility of getting 
trapped into locally optimal points. The accuracy of 
these approaches depended too much upon both 
initially starting points and the degree of system's 
nonlinearity. 

It is better to hybrid these two different approaches 
such that the explorative and exploitive abilities are 
increased simultaneously. Memetic algorithms (MAs) 
are the hybrid programs which incorporate local- 
search methods into global-search technologies. MAs 
are successful in solving various optimization 
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problems of other fields. Harman and McMinn (2010) 
empirically studied the performance of local and 
global search in MAs. Soh et al. (2010) proposed 
archiving molecular and basin-hopping algorithms to 
identify low-energy pure- water isomers. Ahn et al. 
(2010) designed an electromagnetic system. Meuth et 
al. (2009) made a review of Mas and proposed a high- 
order meta-learning system. Kramer (2010) proposed a 
Powell-ILS strategy. Tsoulos and Lagris (2008) 
hybridized these two in series. Yang and Jat (2011) 
proposed a guided-search scheme. 

Most hybrid techniques used evolutionary algorithms 
to identify a suitable initial start for local-search 
learning. In this paper, we propose an inverse aspect: 
Integrate stochastic operations into the gradient-based 
optimization, instead of incorporating the latter into 
the former. This technology avoids from sticking in 
local minima by adjust the climbing direction and 
some stochastic operations. The proposed technology 
was tested with four biological systems. To further 
exhibit the exploration performance the learning was 
implemented in a wide search space ([0, 100] for rate 
constants and [-100, 100] for kinetic orders) with a bad 
initial start (all parameters were initially set to 80- 
neighborbood). 


Gradient-Based Metaheuristics Optimization 

Based on biochemical system theory, the net influx 
( V t + ) and efflux ( Vf ) of an S-system are approximated 
as power-law functions. The concentration change of a 
metabolite, protein or gene is expressed by Eq. (1) 
(Voit, 2000), 


x, = v; - vr = a t n xf j - p t n x ) , <i) 


7=1 


7=1 


for z=l,2,..., n, where n and m are the numbers of 
dependent and independent variables, respectively, on 
and i Si are rate constants, and gij and hi] are the kinetic 
orders that denote the interaction from X] to Xi. 


In order to construct such a highly dimensional 
nonlinear system, a gradient-based metaheuristics is 
proposed. The computational method starts with the 
hill-climbing optimizer, and solves the stagnation 
phenomena by both differential climbing operations 
and migration synchronous evolution. The former 
operation is to adjust the climbing direction and the 
latter is to widen the searching such that a valid 
escape is ensured. Fig. 1 is the proposed flowchart. 



FIG. 1 THE FLOWCHART OF GRADIENT-BASED 
METAHEURISTICS OPTIMIZATION. 


Times Series Data 

The time-courses data sets for training and testing 
were generated from the S-systems, which were cited 
from the published papers (Tsai and Wang, 2005; 
Kikuchi et al., 2003; Voit and Almeida, 2004), in 
different sets of initial conditions. 

Adjust Climbing Direction 

The climbing direction is differentially adjusted 
toward individuals over the searching space: 

I g = I s +r l *d l *O g -I rl ) + r 2 *d 2 *(I g -I r2 ), (2) 

where d\ and d 2 denote the size of a step (d 1 =d 2 = 1 for 
our systems), r x and r 2 are random factors, and I rX and 
hi are two randomly selected individuals. 

Migration Synchronous Evolution (MSE) 

Real-value-coded genetic algorithm is used for the 
gradualness of the variables. We encode the unknown 
parameters of the S-system as a chromosome, and 
arrange all chromosomes in order according to their 
fitness values. The best individual lb (=Ii) has the 
smallest residual error. 

1) Synchronous Mutation 

To ensure global searching, we introduce a 
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synchronous mutation operation. Instead of 
adopting only one gene (one-point mutation), two 
genes (two-point mutation), or a fixed number of 
genes (mask mutation) for mutation, we let all of 
the genes mutate with a probability that is assigned 
by the designer. Those genes with the mutation 
probability q over the randomly given threshold 
r 2 are qualified to replace their original ones: 

r%> i f r \ >r 2 

x ij = < 

Xfj , otherwise 

where Xy is the yth chromosomes in the zth 



individual, and q , r 2 e [0, 1] are two random 
numbers. The synchronous-mutation operation 
brings the population sufficient diversity, but may 
induce a leak in hill-climbing. Therefore, the elitism 
strategy is introduced to compensate the leak. 
Elitism keeps the best-so-far individual to survive 
for each generation and ensures the best 
characteristic to pass down. 
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FIG. 2 SYNCHRONOUS MUTATION. 


2) Migration 

To achieve a wider search we further introduce the 
migration operator such that the new generated 
population covers the entire search space. The 
migration operation is executed in each iteration. 
After this operation new chromosomes are 
generated. The ;th gene of the chromosome I t is 
changed to 

_ \*bj + r 2 x (*/,min - Hj ) if h > n 
ij 1 Xfy + 1'2 x (xj max - x/y ), otherwise 

Xfoj X j 

where A b = — , x bj - are the y'th gene of the 

x j, max — x j, min 

best chromosome, xj ?max and Xj m - m are the upper 
and lower bound of the yth gene, respectively, and 
q , r 2 <e [0, 1] are two random numbers. 



Artificial Experiments and Discussion 

In order to examine the performance of the proposed 
technology, we tested it with a three-gene cascade 
pathway (Tsai and Wang, 2005), a four-gene genetic- 


branch pathway (Voit and Almeida, 2004), a five-gene 
small-scale genetic network (Kikuchi et al., 2003; 
Hlavacek and Savageau, 1996), and a medium-scale 
genetic network (twenty genes) (Noman and Iba, 2006). 
We initialized all parameters at the neighborhood of 
80, and set the range of the rate constants as [0, 100] 
and that of the kinetic orders as [-100, 100]. All 
computations were performed on an Intel Core duo 
3.16 GHz computer using Microsoft Windows XP. The 
performance is determined by the weighted mean- 
square error (residual error). 


E 


i N 

-v. 

N “ 

i=i 


max(x* ) 

exp J 


(5) 


where x [ and x l , z=l,..., n are, respectively, the 

estimated concentration and artificial-measured 
concentration, t a is a time-weighting factor, and N is 
the number of sampled data points. The cubic-spline 
technology was used to generate the smooth profile of 
the time-series data 


Cascade Pathway System 
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FIG. 3 A CASCADE PATHWAY (TSAI AND WANG, 2005) 

We first consider the cascade pathway shown in Fig. 3 
(Tsai and Wang, 2005). This is a three-step system with 
two feedback signals. The system has three dependent 
variables ( x x , x 2 and x 3 ) and one independent 
variable ( x 4 ). The respective S-system is 

xi = 6ir 1 x| 12 xf 13 x 4 - P x x \ n , 

xi - a 2 xf 21 - P 1 x h 1 11 , (6) 

X 3 = Cl^xf 32 - P 3X3 33 . 

The values of parameters are listed in Row 'True" of 
Table 1 We generated eight-set artificially 
experimental concentration data for the training. The 
simulation time is 8 seconds, and the sampling time is 
0.02. After the training, we get the estimated 
parameters, which are shown in Row "simulation" of 
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Table 1. The estimated residual error in Eq. (5) is 
1.0591741E-10. 

Genetic-Branch Pathway System 

The second system is the genetic-branch pathway in 
Fig. 4, which is used by Voit and Almeida (2004). The 
respective S-system is described as Eq. (7). 

xi = a x x 3 u x 0 - fi x x \ n , 

x 2 =a 2 xf 21 ~P 2 xl 12 , 

X3 =a 3 X 2 32 -P 3 X 3 33 X^ 34 , 

X4 =a 4 xf 41 - p 4 x\ 44 . 

The system has four dependent variables (x 1 ,x 2 ,X 3 
and x 4 ) and one independent variable ( x 0 ). Row 
"True" of Table 2 shows the parameters. Eight sets of 

TABLE 1 TRUE AND ESTIMATED PARAMETERS OF AN S-TYPE SYSTEM FOR A CASCADE PATHWAY (3 GENES) IN FIG. 3. ROW "TRUE" IS THE PARAMETERS OF A 
TRUE S-SYSTEM. ROW "SIMULATION" IS THE ESTIMATED PARAMETERS FOR A WIDE SEARCH SPACE ([0,100] FOR RATE CONSTANTS AND [-100, 100] FOR 

KINETIC ORDERS) WITH A BAD INITIAL START (80 FOR ALL PARAMETERS). 
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XI 
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-0.1 -0.05 0.5 

true 

X2 
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1.44 

0.5 
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X3 

3 

7.2 


0.5 0.5 


XI 

9.9999645E+00 4.9999930E+00 


-9.9999021E-02 -5.0001491E-02 5.0000005E-01 

simulation 

X2 

2.0000463E+00 1.4400463E+00 4.9999122E-01 

4.9998848E-01 


X3 

3.0002409E+00 7.2001110E+00 


4.9996317E-01 4.9996313E-01 


TABLE 2 TRUE AND ESTIMATED PARAMETERS OF AN S-TYPE SYSTEM FOR A GENETIC BRANCH PATHWAY (4 GENES) IN FIG. 4. ROW "TRUE" IS THE 
PARAMETERS OF A TRUE S-SYSTEM. ROW "SIMULATION" IS THE ESTIMATED PARAMETERS FOR A WIDE SEARCH SPACE ([0,100] FOR RATE CONSTANTS AND 
[-100, 100] FOR KINETIC ORDERS) WITH A BAD INITIAL START (80 FOR ALL PARAMETERS). 



artificial experimental data ( x l exp ,i = 1 are used for 

the training. The simulation time for each experiment 
is 8 seconds and the sample time is 0.02. Row 
"Simulation" of Table 2 shows the estimated 
parameters, all of which are nearly the same as their 
respective true values. The mean-square error in Eq. (5) 
is 2.2096583E-11. 
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Small-Scale Genetic Network 

We further consider a small-scale genetic network, as 
shown in Fig. 5 (Kikuchi et al., 2003; Hlavacek and 
Savageau, 1996). The diagram shows that the 
transcription of the gene x 6 is regulated by two 
feedback signals from x 3 and x 4 , respectively. The 
following is the respective S-system: 

xi = a x x f 13 xf 15 x 6 - p x x x n , 

X 2 = a 2 xf 21 x 7 - P 2 X 2 22 J 

X 3 = a 3 xf 32 x 8 - P 3 x 2 32 x 2 3 \ (8) 

V . _ v S 43 v &45 v P y.^44 

JL4 — LC 4^3 Ji 5 Ji ft JJ 4 A 4 , 

X 5 = ^ 2*^4 54 W P 7 ^5 ^ • 



FIG. 5 A SMALL-SCALE GENETIC NETWORK (KIKUCHI ET AL., 
2003; HLAVACEK AND SAVAGEAU, 1996) 

There are five dependent variables x it i = 1 ,..., 5 and 
three independent variables x t , i= 6, 7, 8. The values of 
the rate constants and kinetic orders are listed in Row 
'True" of Table 3. We generated the same numbers of 
data sets as the first and second systems. Eight-set 
artificially experimental data X ^ xp = were 

generated from the S-type dynamic system. The same 
technology was used to get the smoothing profiles of 
these eight-set data sets. Artificial experiments were 
proceeded from the time-instant t=0 to t=0.5 sec. with a 
sampling time 0.0125. The estimated parameters are 
shown in Row "Simulation" of Table 3. The mean- 
squared error in Eq. (5) is 9.3606730E-10. 

Medium-Scale Genetic Network 

Finally, we consider a twenty-dimensional system, as 
shown in Fig. 6. The medium-scale genetic network 
was used by Noman and Iba (2006). The system has 
twenty dependent variables x if 20, but not an 


independent variables exists in the system. The 
degradation rate (the efflux V[ ) of each constitute ( x t ) 
depends only on himself, which is drawn as the self- 
feedback signal, as shown in Fig. 6). The regular 
signals are denoted by the directed branches. From the 
diagram, we know that the production rate of the 
constitute x 7 depends onx 2 , x 3 andx 10 . The respective 


S-system is 
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FIG. 6 A MEDIUM-SCALE ARTIFICIAL GENETIC NETWORK 
(NOMAN AND IBA, 2006) 

The rate constants and kinetic orders are listed in 
Column "True" of Table 4. Eight sets of artificially 
experimental data (xg Xp ,/ = l were generated. The 

simulation time is 1.8 seconds and the sampling time 
is 0.01. The estimated parameters are shown in 
Column "Simulation" of Table 3. The mean-squared 
error in Eq. (5) is 3.5860070E-08. 
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TABLE 3 TRUE AND ESTIMATED PARAMETERS OF AN S-TYPE SYSTEM FOR A SMALL-SCALE GENETIC NETWORK (5 GENES) IN FIG. 
5. ROW "TRUE" IS THE PARAMETERS OF A TRUE S-SYSTEM. ROW "SIMULATION" IS THE ESTIMATED PARAMETERS FOR A WIDE 
SEARCH SPACE ([0,100] FOR RATE CONSTANTS AND [-100, 100] FOR KINETIC ORDERS) WITH A BAD INITIAL START (80 FOR ALL 

PARAMETERS). 
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TABLE 4 TRUE AND ESTIMATED PARAMETERS OF AN S-TYPE SYSTEM FOR A MEDIUM-SCALE GENETIC NETWORK (20 GENES) IN FIG. 6. "TRUE" IS THE 
PARAMETERS OF A TRUE S-SYSTEM. "SIMULATION" IS THE ESTIMATED PARAMETERS FOR A WIDE SEARCH SPACE ([0,100] FOR RATE CONSTANTS AND [-100, 

100] FOR KINETIC ORDERS) WITH A BAD INITIAL START (80 FOR ALL PARAMETERS). 
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Discussion 


Learning in the wide range and the bad initial start is 
to show that the exploratory ability of the gradient- 
based method is largely increased. We observed that 
the proposed technology is able to escape from local 
minima and keeps the advantage of computation time. 
The results of Tables 1-4 show the high accuracy of the 
proposed method. The mean-square error in Eq. (5) is 
1.0591741E-10 for 3-gene, 2.2096583E-11 for 4-gene, 
9.3606730E-10 for 5-gene, and 3.5860070E-08 for 20-gen 
systems. 

We further use Figs. 7-10 to show the advantage of 
convergence. We observe that most of the curves 
converge at around 1000 fitness evaluation. In 
addition to the strongly nonlinearity, biological 
systems are always stochastic wherein the time-series 
data are heavy noisy. Therefore, we further consider 
the artificial data with 10% random noise. Fig. 11 
shows that the proposed method is able to predict the 
dynamic behaviour of the system even in a noisy 
environment. The results are comparable to Fig. 4 of 
Sriyudthask's paper (2013) where 0 to 5% noisy data 
are used. 
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FIG. 7 THE CONVERGENCE OF THE PROPOSED METHOD 
APPLIED TO THE CASCADE-PATHWAY SYSTEM IN EQ. (6). 
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FIG. 8 THE CONVERGENCE OF THE PROPOSED METHOD 
APPLIED TO THE GENETIC-BRANCH PATHWAY IN EQ. (7). 
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FIG. 9 THE CONVERGENCE OF THE PROPOSED METHOD 
APPLIED TO THE SMALL-SCALE GENETIC NETWORK IN EQ. 

(8). 
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FIG. 10 THE CONVERGENCE OF THE PROPOSED METHOD APPLIED TO THE MEDIUM-SCALE ARTIFICIAL GENETIC NETWORK IN 


EQ. (9). 

cascade pathway with 1 0% random noise 



FIG. 11 ROBUST EXAMINATION FOR THE GRADIENT-BASED METAHEURISTICS OPTIMIZATION IN THE CASCADE-PATHWAY. DOT 
POINTS ARE DATA WITH 10% RANDOM-NOISE CONTAMINATE. SOLID CURVES ARE THE ESTIMATED PROFILES. A WIDE SEARCH 
SPACE ([0, 100] FOR RATE CONSTANTS AND [-100, 100] FOR KINETIC ORDERS) AND A BAD INITIAL START (80-NEIGHBORHOOD FOR 
ALL PARAMETERS) ARE USED. INITIAL CONDITIONS ARE 20% BEYOND THE TRAINING RANGE. 
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Conclusions 

Identifying a dynamic biological system from time- 
series data is a central theme in systems biology. S- 
systems are demonstrated to be a good approximation 
of smooth nonlinear systems. However, the parameter 
identification of S-systems is challenging because the 
S-system is described as coupled highly nonlinear 
differential equations. How to make a trade-off 
between the accuracy (reliability) and computation 
time is important. Instead of improving the local- 
search ability of the population-based computational 
methods, we propose an inverse aspect: The 

incorporation of stochastic-search operations 
(migration synchronous evolution) into traditional 
gradient-based optimizers. The synchronous mutation 
is to increase the population diversity and the 
migration operation to widen the searching. The 
simulation results exhibit that the proposed scheme is 
able to escape from local minima in a reasonable 
computation cost, even the dimension of the system is 
as high as 20 and the searching space is wide. 
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